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Pl, Abstract 

We discuss the interpolation of the electric and magnetic fields within a charge- 
conserving Particle-In-Cell scheme. The choice of the interpolation procedure for 
the fields acting on a particle can be constrained by analyzing conservation of 
the energy and the particle generalized momentum. The better conservative prop- 
erties are achieved, if the alternating-order form-factor is used for interpolation, 
i/i • which combines the lower-order and higher-order interpolation from integer and 

semi-integer points of a staggered grid. This approach allows us to significantly 
0^ improve both the results quality and the computational efficiency for the charge 

conserving scheme. 

Keywords: Particle-In-Cell, conservative scheme, charge-conserving scheme 

1. Introduction: charge-conserving PIC schemes. 

Here we discuss the conservative properties of the Particle-In-Cell (PIC) nu- 
merical schemes. Recall that usually the conservative schemes are employed to 
solve the system of conservation laws. For example, the continuity equation, 
+ V • (pmu) = 0, for the fluid mass density, p^, may be advanced through the 
time step. At, using the conservative scheme: Vi{pm)'i'^^ = Vi{pm)i — At J2j cnj, 
with the computational domain split into a set of control volumes ("cells"), the 
mass density at a given time instant, = nAt, averaged over the volume of the 
cell, i, and the mass flux through the ij face averaged over the time step: 

If 1 f 

{pm)7 = yj iPm)t=t"dVi, aij = — J^^ dt J {dSij ■ up„), (1) 

where u is the fluid velocity and the face area vector, Sij is directed from cell i, to 



cell j. Since cr,y = —o-ji, the total mass is conserved: J2i iPm)i = J2i (p 



mJi 
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This idea may be employed to assure the charge conservation within the PIC 
method, although the governing equations differ from conservation laws: 
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where the index p enumerates particles (electrons, ions), Wp is a momentum-to- 
mass ratio, and other notations are usual. The particle charge density, p, obeys not 
only the continuity equation, || + V ■ J = 0, but also the Poisson equation. There- 
fore, the charge conservation property is formulated as the relationship between 
the charge density and the electric field, E, and/or the electric current, J. 

In the present paper we discuss only the CHarge-Conserving PIC (ChCPIC) 
schemes. The need to introduce this family of numerical schemes is that to simu- 
late fully relativistic motions of plasmas in strong rapidly varying electromagnetic 
fields (such as the laser field, for example), one may want to directly advance the 
electric field using the first Eqs.©. The "charge conservation" property in this 
case allows one to fulfill the Poisson equation for this advanced field without 
solving the equation directly. The way to assure that the charge conservation was 
developed in IE 0, S, H, Hi (see alsojp, 0, S, @] and references therein). First, a 
staggered grid should be used as in [|lO|] to ensure the finite-difference approxima- 
tions for [V X B], [V X E] to be divergence-free, as we discuss briefly in Sec. 2. 
Second, the currents through the cell faces should be calculated in such way that 
their divergence balances the charge leakage from the cell. In Sec. 3, this is done 
using the virtual path integration to solve the time integral in Eq.©. However, 
once the way to compute the particle current is modified, we should also modify 
accordingly the scheme for interpolating electric and magnetic fields acting on 
this particle, which is the goal of this paper. 

Note an important distinction of the ChCPIC schemes from the general cloud- 
in-cell scheme, which is explicitly pronounced once we follow the conservative 
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schemes ideology. On one hand, a "cloud" within the framework of the cloud- 
in-cell scheme can be thought of as some charge density distribution, pp(x, Xp), 
centered about the particle coordinate vector, Xp. Following this ideology, the 
point-wise values of the shape-function (form-factor), pp(x, Xp), at the points x, 
where the electromagnetic fields are localized, should be used then as the inter- 
polation coefficients to average the electromagnetic fields acting on the particle, 
and to calculate the Lorentz force. Contrarily, the conservative schemes approach 
assumes that the cloud charge density should be integrated over some control 
volume or its faces, to represent the contribution to the plasma density and the 
plasma electric current (see [U), thus making inapplicable the form-factor point- 
wise values. To constrain the field interpolation procedure, in Sec. 4 we discuss 
the accuracy of conservation, for the energy integral and for the particle general- 
ized momentum. Their conservation in the ChCPIC scheme can be achieved, if 
the alternating-order form-factor is employed in the interpolation procedure. This 
alternating-order form-factor combines the interpolation of different order for dif- 
ferent physical variables to interpolate. Equivalently, it integrates the "cloud" 
charge density over faces and edges of the control volume to interpolate face- 
centered and edge-centered electromagnetic field values. 



2. Grid geometry and notations. 

We use a 3D Cartesian grid in the domain < x < L^, < y < Ly, 
< z < Lz, split for * Ny * cells. The coordinates of the cell corners are 
(i Ax, j Ay, k Az), where i, j, k are integers and Ax = L^/N^, Ay = Ly/Ny, 
Az = Lz/Nz are the cell sizes. Then we introduce the normalized coordinates, 
X = x/Ax, y = y/ Ay, z = z/ Az, and time, i = t/ At, and use them below 
with no tilde. Magnetic field, electric current, and particle momenta are defined 
at semi-integer time instants, t = n + 1/2, the electric field and the particle co- 
ordinates - at integer time instants, t = n. In the normalized coordinates, Eqs.© 
read: ^ 

= — ■ diag{cr,,Cy,Cz), p = ^|7(5(x-Xp), (7) 
dt c p ^ 

where V = AxAyAz and = cAt/Ax, .... The grid functions are defined: the 

(x) 

cell-centered charge density, pi+i/2,j+i/2,k+i/2; the electric field, -E'-/+i/2,fc+i/2' 

(v) (z) (x) 

^i+i/2j,fc+i/2' ^i+i/2,i+i/2,fc' and the current density components, J> 

j+l/2,fc+l/2' 

"^1+1/2 j fe+1/2' "^1+1/2 j+1/2 k' defined at the centers of the faces, normal to the axis, 
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X, y, z\ and the magnetic field components, 5|^i/2j- fe, ^ljVi/2,fc' ^ij,k+i/2' 
fined at the midpoints of the edges directed along the axis, x, y, z. The subscript 
indexes denote coordinates of the point at which the grid function is defined. 

The PIC scheme as taken from [lllll with the suggested modifications is de- 
scribed in the Appendix. The algorithm involves interpolation for the fields acting 
on a particle: 
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' n\ 


^i,j+l/2,fc' 




(x;) = 


E f^iJ,k+l/2 
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where a, (3 are the weights, their sums for a given particle should be equal to one: 

E «SViA.+i/2(xp) = 1, ... E P\%2,M) = 1' - (8) 

Herewith, we provide only the expressions for x component of vectors, whenever 
possible, denoting the generalization for the other components by The con- 
tribution to the current density from a charged particle can be expressed in terms 
of the particle position, x^, and its velocity, u^+^Z^, and/or Xp+^: 

(x) n+l/2 _'<^(lp^{x) , n ^n+l\ 



j(x) n+L/2 _ ypp 
'^i,j+l/2,k+l/2 ~ T/?; 



j+l/2,fc+l/2 - T7Sjj+l/2,fc+l/2l^p 



x„, x. 



p 



j{y) n+l/2 _ (IpAy) /^n „n+lN 

"^i+l/2j,fc+l/2 ~ T/'^j+l/2j,fc+l/2V^p' ^p 



P 



V 



j{z) n+l/2 _ sr^ QpM) / n „n+lN /qn 

•^j+l/2j+l/2,/c ~ T^'^i+l/2j+l/2,fcl-^p' /• W 

P ^ 
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The advantage of the staggered grid is that the magnetic field divergence does 
not change and equals zero as long as it is initially equal to zero. Analogously, 
[V X B] term does not affect the electric field divergence. The Poisson equation, 

47rpr+i/2,7+i/2,fc+i/2 = (-^i+i,7-+i/2,fc+i/2 ^ ^i}+i/2,k+ii2)l + Satisfied, if: 



i+l/2,j+l/2,fc+l/2 ~ l-^i+l,i+l/2,fc+l/2 ^i,j+l/2,k+l/2) 

_ j{x) n+1/2 

) ^'^1 /9 J- -1-1 /9 ^ 



.„ _ „+l jix) n+1/2 _ j{x) n+1/2 

ri+l/2,j+l/2,k+l/2 Pj+l/2j+l/2,fc+l/2 _ j+l/2,fc+l/2 i,j+l/2,k+l/2 



At Ax 

j{y) n+1/2 _ jiy) n+1/2 Jz) n+1/2 _ j{z) n+1/2 

"'i+l/2,j+l,fc+l/2 "'»+l/2,i,fc+l/2 '^i+l/2j+l/2,fc+l '^i+l/2,j+l/2,k 

Ay Ax • ^ ^ 



3. Charge density and charge conservation law. 

3.1. Form-factors. 

To discretize the charge and current densities, one needs to specify the nu- 
merical representation for S functions in Eqs.©. We do this using the family 
of form-factor functions, f^^\x,Xp), where the form-factor of a zero order is a 
cap-function: f^^\x,Xp) = 1, and the higher-order form-factors are recursively 
defined: f^''^^\x,Xp) = /^J^jyg^ /^'^(x', All form-factors: (1) are symmet- 

ric functions of x — Xp-, (2) turn to zero at \x — Xp\ > (/ + l)/2; and (3): 

. yj!^ . /.'.(. + 1/2. - - 1/2. 

We are interested both in point values of the form-factor function, and in its 
integrals over the grid size. So, for a chosen form-factor, /(x, Xp) = f^''\x, Xp), 
we introduce: 

fiip^p) fih-^p)i Fiip^p) / /(-^ , Xp)dx 

J —oo 

and 

ri+l 

AFj+i/2(xp) = Fi+i{xp) - Fi{xp) = J f{x', Xp)dx'. (11) 

By definition, AFi^ii2{xp) = f^^~^^\i + 1/2, Xp). The applicability of the form- 
factors for constructing the interpolation weights, which should satisfy Eq.®, is 
ensured by the equation: 

J2f'\x + t,Xp)= I f-'\x,Xp)dx = l. (12) 

~ J\x-Xj,\<l/2 
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From Eq.(fT2l) we can obtain yet another identity to be used below: 

=x;-a:p. (13) 

i 

To prove Eq.(fT3]) one can note that for Xp = x'^ both sides of the equation turn to 
zero and that on taking the derivative of this equation over x'^ we obtain the earlier 
proven Eq.(fT2l). 

3.2. Conservative scheme for electric charge 

A particle can be thought of as a cloud with the charge density, Pp(x, y, z) = 
Qpfix, Xp)f{y, yp)f{z, Zp) /V . Following the conservative scheme idea, we define 
the contribution from this particle to the charge density grid function not as a point 
value of Pp in the cell centers, but via Eqs. dlllll) : 

Pr+i/2,,+i/2,;t+i/2 = ^ E gpAF,+i/2(x^)AF,H.i/2«)AFfc+i/2(Zp"0. (14) 

Note that the integrated over the cell size form-factor value, AFj+i/2(a;p), ... is at 
the same time the point-wise value of the form-factor function of by unity higher 

order, AF,^.l/2(x^) = /^^t/2(^p)■ 

Assuming that within the time interval, (n, n + 1), the particle moves from 

the point to the point Xp"*"^ along an arbitrary virtual path x„(t), we define the 

particle currents, ^(^'J^'^), in Eq.® following Eqs. dlll 11141) ): 



n+l i+1 fe+1 



d-ViA^li/2(x,",x,"+^)= / / J ^^f{x^)f{y\y.)fiz',zMz'dy'dt 



n+l 



Ax J dFiixy] 
'At J 



^^AFj+i/2iy.)AFk+i/2iz„)dt, (15) 

n 

dSS£i/2(Xp,xr^) = 7 ^^AF.+,/2(x.)AF,^,/2(.„)cit, (16) 

n 

ea/t>WK'-r') = -|^ / ^^AF.+,/2(x.)AF,^,/2(y.)rft. (17) 

n 

By definition, Eqs. dMHfTI) satisfy the charge conservation law as in Eq.(fTOl). This 
can be verified by observing that the linear combination of Eqs. dlSHFTI) as pre- 
sented in Eq.([IOl), reduces to / ^[AFj+i/2(x^) AFJ+l/2(^/^,)AFfc+l/2(^^)]rf^. 
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Discuss a choice of the virtual path. For the straight path, x^(t) = + 
{t — n)Up+^/^ ■ diag{cx, Cy, Cz) /c, the integrands in Eqs. (|15H17l) are piecewise 
polynomials of the order of 3/ + 2. The integration was performed in |0],[l5l] 
only for the lowest order form-factor, / = 0. With a finite time step, especially 
for higher-order form-factors, the integration becomes sophisticated, because the 
form-factor is a piecewise smooth function. All ranges of smoothness should be 
accounted for separately while evaluating the integral analytically. 

To overcome this difficulty, we observe that the integrals in Eqs. (|15H FT] ) can 
be solved, if the virtual path is composed of the edges of the rectangular box, such 
that the points :>i.^\ x^,"^^^ are the opposite corners of this box and its edges are 
parallel to the coordinate axis. Upon calculating the integral in Eq. (fT5l) as a sixth 
of a sum of the integrals along six possible virtual paths, we find: 



Ax) n+l/2 , n n+U _ 
^i,j+l/2,k+l/2y^p^^p )- g 



Fii^p ) ^ii^p] 



+ 

1 Ax 



X 
i'<i 



x{2 [AF,+i/2«+^)AF,+i/2(^r') + AF,+v2«)AF,+i/2(;^; 
+AF,+v2(|/r')AF^+i/2(^;) + AF,H.i/2(|/,")AF,+i/2(z;+^)} = 

E AFi:l/,(x,)[AF,(:;/,(y,)AFi+i/,(.,) + \AF^;Uyp)^Hl\/2i^p)l 
i'<i 

Ay) n+l/2 , n n+l^ _ V- A p(-) /„, ^^ 

Q+i/2,i,fc+i/2l^p' )- A A. 2^ ^^j'-i/2yyp)^ 

j'<j 

Az)n+\I2 / n n+lN _ V- A i^(-) / N 

?i+l/2j+l/2,fcl-^p'-^p )— AX. ^-^fc'-l/2l^pJ^ 

fc'<fc 

X [AF,^)/2(x,)AF,(+;/2(y,) + iAF/-)/2(x,)AF,y/2(y,)]. (18) 

Here AF^JJ/^I^/p) = AF,+i/2(y;+i) ± AF,+i/2(i/p, .... A recursive formula, 
Fi{xl'^^)-F,{xl) = F,_i{x'^+^)-Fi^i{x;) + AFl:l^{xp), allows us to calculate 

Fiix"^^^) - Fiix'l) = J2i'<iAF^r\^^{xp). The scheme as in Eq.® is not new 
and was obtained from different considerations and in a different form in [6] (see 
also discussion in [[siS])- 
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Note a useful property of the interpolation coefficients. Summing up the con- 
tributions to the electric current at different faces and using Eqs. (ll2ll3l ) we obtain 
the following exact relationship: 

V- ^(x) n+1/2 / n n+U _ (x) n+1/2 
Sj,j+l/2,fc+l/2l^p' ) ~ "'p > 

and the analogous relationships for y and z components. These relationships for 
the interpolation coefficients for the electric current are similar to those for the 
electric and magnetic fields (see [8]). 

4. Interpolation procedure for the electric and magnetic fields. 

To interpolate the charge density we had to integrate the form-factor over the 
cell volume. Here we show that to interpolate the electric and magnetic fields the 
form-factor should be integrated over the cell faces and edges correspondingly. 

4.1. Energy integral 

Consider the energy, £p, of the system consisting of the electric field and a sin- 
gle charged particle. Both the particle energy and the electric current contributing 
to the electric field energy change are additive by particles, hence, so is the error 
in the energy conservation and thus, we can calculate it for a single particle. For 
simplicity, we assume zero magnetic field because this field does not affect the 
particle energy and any change in the magnetic field energy is balanced with that 
for the electric field. At time instant, t = n, the energy, S^, equals: 

W(wr^/^)^ + c2 + i^ (^tWn-i/^)' + ...] +^(ur^/^-E(x;)), 

the last term is to advance the particle energy through a half time step. The change 
in the particle energy with the use of Eq.(|33l). can be approximated as: 

v/(wr^/^)2 + c2 - v/(wr ^/^)2 + c2 ^ 1^ (e(x;) . {u;-y^ + u^^/^)) , 

with the error, 0((At)^). Neglecting this error, we derive the change in the total 
energy, using Eqs. (l9l35l) : 

^r' - ^; = f f E Mi^:'r - (41:)^] + ■] + 
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^((E(x^) + E(x^^)).<+V2 



QpAt 



(E(x- 



\^ ( A c{^}n+lf2 . p{a;)n+l/2 , 



where: 



(x) n+1/2 _ ^ ^{x) n+1/2 
jj + l/2,fc+l/2 ~ ^I'^x face 



qpA.t , „(x) n 



u 



(x) n+1/2 



X face 



(x) 



X faceV-'^p 



4 

n+l' 



face 



-^x face J ^ 



2^ (x) n+1/2' 
"^x face 



. n+1/2 
^3^x face 



7T {qpMf' ^(a;) n+1/2 



V 



(x) n+1/2 (Ax) 



-C 

Sx 



face 



X 



face / face 



n+1^ 



"x face" stands for i, j + 1/2, A; + 1/2. If the interpolation weights for the electric 
field match those for the current in the way as follows: 



(x) n+1/2 r (a;) 



[«rLe(x;) + «?L(x;+^)] = 2^i^ 



■■{x) n+1/2/ n n+l 



face 



Xp, Xp 



(19) 



then Ai£ = and the energy defect is small: A^S ~ (At)^. 

The problem is that, in contrast with an ordinary PIC scheme, within the 
ChCPIC scheme one can hardly satisfy (fT9l ) exactly. Eq.(fT9l) can be obtained 
as the trapezoidal estimate for the integrals in Eqs. (|15H17] ). if the interpolation 
weights, a, are chosen as follows: 



"j,j+l/2,fc+l/2V^p> 

(T^^^ fx"! 
"i+l/2j,fc+l/2l^p^ 

"j+l/2j+l/2,fcl^p> 



AF,+i/2(x;)/,(y;)AFfcH_i/2(^;), 



(20) 



The accuracy of Eq.(fT9l) and the energy defect while using Eqs. (1 15m 8 1201) are 
controlled by the choice of the form-factor order. Using the estimate: 



At /\t r^~x'^ 

— [i^K^p) - ^.«')] = ^ /(^ - xp)d{x - 



Ax r~^v 
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^{x) n+1/2 



1 



12 

where we substituted 



At /-i/a df{x-Xp) 



Ax J-i/2 dx 



x-Xp = t-{x;+x;+^)/2+g{x;+^-x;) = t-{x;+x;+^)/2+gul'''> ^+'/^At/Ax, 

and noting that for / = 0, 1,2 correspondingly the form-factor, f{x — Xp), its 
first or second derivative are bounded, we can evaluate the accuracy of the en- 
ergy conservation using the above estimate: Ai£ ~ 0((At)'"''^). If / = or, 
alternatively, if Eq.(fT9l) is not fulfilled at all, the energy does not conserve in the 
ChCPIC scheme: AiS ~ 0{At). We conclude that both the use of higher-order 
form-factor (/ > 1) and the interpolation following Eq.(l20l) are desirable. 



4.2. Generalized momentum conservation 

The governing equations (I!® conserve the projection of the generalized par- 
ticle momentum, m.pWp + qpA/c, on the given direction g, if the electromagnetic 



field is constant along this direction, i.e. (g ■ V) A = (see lll2n ). Indeed, 



^(nipWp + ^A) = ^ ([u, X [V X A]] + (u, ■ V)A) = ^ 



and g • -^{mpWp + qpA/c) = as long as (g ■ V)A = 0. To verify the gen- 
eralized momentum conservation within the ChCPIC scheme, the latter should 

(x) 

be formulated in terms of the vector potential. The grid functions ^^^4.1/2 fc+1/2' 
4+i/2,j,k+i/2' ^m/2,i+i/2,fc are introduced at the same points as 4?+i/2,fc+i/2' 

^i+1/2 j k+i/2' ^i+1/2 j+1/2 k ^^'^ '•h^ ^^^^ derivative of the vector potential grid 
function may be expressed in terms of the electric field: 

A (x) n+1/2 _ .(x) n~l/2 _^A + p(2:)"' 
^ij+l/2,fc+l/2 ^ij+l/2,fc+l/2 '"^''-^ij+1/2, fc+1/2' 

.(jy) n+1/2 _ .(y)n-l/2 _„A:^P^(^)" 
^i+l/2j,fc+l/2 ~ ^j+l/2,j,fc+l/2 ^^''-'^i+l/2j,fc+l/2' 

.(2) n+1/2 _ .(2) n-1/2 _^Ay-p(^)" nO^ 
^i+l/2j+l/2,fc ~ ^i+l/2j+l/2,fc ''^''-^i+l/2,j+l/2,fc' V^^^* 

in accordance with the first of Eqs.®. The electric field and the vector potential 
at X = Xp are also linked via this equation, hence, the vector potential in a particle 
location, A{t, Xp) should be interpolated with the same weights, a, as we use for 
interpolating the electric field: 

y4^^^(t, Xp) = X! A-j+l/2,fc+l/2'^!j+l/2,fc+l/2(^p)' ••• (^3) 
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On the other hand, the magnetic field acting on the particle can be interpo- 
lated via the grid function of the magnetic field, B^"^^ = J2i,j,k Pi+1/2 j k^i+1/2 j k- 
In turn, the latter grid function can be expressed in terms of that for the vector 
potential, as the discretization of the second of Eqs.®, on the staggered grid: 

^i+l/2,j,k \'^i+l/2,j+l/2,k ^i+l/2j-l/2,fc 



1 



Ay) _ Ay) 

■r^A I 1 /o „• I. I 1 /o 



SO that: 



JD [^L, }^) — \yi+l/2,j,k yi+l/2,j+l,k)^i+l/2,j+l/2,k 

1 



if^i+l/2,j,k /^i+l/2,i,fc+l)A^+l/2,i,fc+l/2' ••• (24) 

Now, we introduce the generalized momentum, 

P-V2 = rripw;-'/^ + ^A{n - 1/2, - w"-^/^— ). 

Evaluating the difference P"+i/2 — ^j^^^ ^j^g transformation as in 

Eq.(l2TI). which allows the generalized momentum conservation, is possible, if the 
differential equation, B = V x A, ?5 exactly fulfilled with the interpolated values 
of the magnetic field, and vector potential, i.e. 

B(xp) = [VpX A(t,xp)], (25) 

where Vp = (-^^f-, ••)• Applying the operator, VpX to Eq.(l23l) and comparing 
the result with Eq.fe4l) we find that Eq.(l25l) is fulfilled if the following set of 
equations holds: 

oix) _ n{x) _ t>'«i+l/2,j+l/2,fc 

f-'i+l/2,j,k f-'i+l/2J+l,k Qy^ ' 

n{x) _ n{x) _ ^"■i+l/2,j,fc+l/2 

^^i+l/2j,fc f^i+l/2,j,k + l ' 



P 



etc. These equations as well as Eq. (|20l) dictate the following choice for the mag- 
netic field weights: 

/5j?i/2,,> = AF,+i/2K)/,(|/p")/.(z;), = /.(xpAF,+v2(2/;)/.(^;), 

/3Sfc+i/2 = f^{xl)f,{yl)^F^^r|2{zl). (26) 
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4.3. Electrostatic limit and a link to the finite element approach 

To discuss the relationship with the important concepts of the momentum con- 
servation and self-force (see detail in[l 1]) we now consider the electrostatic limit 
in which the fields acting on a given macroparticle can be assumed steady-state. 
Within this approximation, the electric field may be expressed in terms of the 
scalar electric potential, $: 

E = - V$ (27) 

The grid function for the scalar potential is cell-centered, $i+i/2,j+i/2,A:+i/2- The 
relationships for the face-centered electric field components are as follows: 

-^ij+l/2,fc+l/2 = \^i-l/2,j+l/2,k+l/2 - $i+l/2,i+l/2,fc+l/2 j , ••• (28) 

The Poisson equation should be fulfilled in the following form: 

— $j-l/2j+l/2,fc+l/2 + 2$j+i/2j+l/2,A:+l/2 " $i+3/2j+l/2,fc+l/2 _ 

(A^P ^ ■■■ ~ 

= 47rpi+i/2j+i/2,fc+i/2- (29) 

With this observations, the contribution from the electric field to the energy inte- 
gral may be transformed using Eqs. (|28|29l) as follows: 



^ I 51 (-^j,i+l/2,fc+l/2)^ + ••• I ~ 2" ^ Pi+l/2j+l/2,fc+l/2$i+l/2,i+l/2,fe+l/2- 
i,j,k I i,j,k 



The key point is that in the last sum the charge density can be expressed in terms 
of the contributions from the given particles, then resulting in the field energy 
transformation to the sum over all charged particles: 

17 XI P»+l/2j+l/2,fc+l/2$i+l/2j+l/2,fc+l/2 = - ^qp^i^ip), 
^ i,j,k ^ P 

where the static potential at the particle location is introduced as follows: 

^i^p) = ^^i+l/2(a;p)^^i+l/2(2/p)AFfc+i/2(Zp)$i+l/2,j+l/2,fc+l/2- (30) 

If the energy is conserved, then the energy integral may play the role of the Hamil- 
tonian function, so that the electric force acting on a particle can be expressed in 
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terms of partial derivatives of the Hamiltonian function over particle coordinates. 
As the result, the electric field acting on the particle should obey Eq.dTTI): 



E(xp) = -Vp$(xp). 

But this is achieved with the choice of Eq.(|20l) for the weight coefficients, to 
interpolate the electric field acting on a particle! Indeed, 

1 d 

~'Kxdx~ ^ ^-^*+l/2(^p)^-^i+l/2(Z/p)AFfc+i/2(2p)$i+l/2,j+l/2,fc+l/2 = 

= ifi+d^p) - fi{Xp))^Fj+i/2iyp)AFk+l/2{Zp)^i+l/2,j+l/2,k+l/2 = 

"ij+l/2,fc+l/2-'^ij+l/2,fc+l/2 ~ -^p K-'^PJy ■■■ 

Here, the conclusion of Langdon is reproduced, that the "energy-conserving" 
plasma simulations models require the contribution from the potential electric 



field to the force acting on a particle to be a gradient of the grid potential (see 11311 '). 
From here yet another way to arrive at the suggested scheme for interpolating the 
electric and magnetic force may be recalled, specifically, the finite elements. The 
nodal amplitudes for fields and the form-factors are equivalent to local polynomial 
representation for the field using finite elements. The introduction of incomplete 
polynomial face (div-conforming) and edge (curl-conforming) elements is gener- 



ally attributed to Nedele c [|14l1 . For the "energy-conserving" scheme presented by 
Lewis in {]3] (see also [16]) the use of such finite elements allows maintaining 
some of the Maxwell equations as the exact relationships between the finite ele- 
ments for the fields. In the application to arbitrary curvilinear grids made in ^ 
this approach resulted in the same form-factors as in Eqs. (l20l26l) . the form- factor 
order being 1 = 1. 

Nevertheless, the new derivation of Eqs. (l20l26l) for newer ChCPIC deserves 
the attention, moreover that the formulation of the conservation laws for such 



schemes is different from the formulation of 11511 . We also emphasized that for 
the form-factor order of / = 1 there is no point in applying the interpolation 
scheme of Eqs. (l20l26l) . as long as the energy (nor the generalized momentum) is 
not conserved anyway with this form-factor. 

4.4. ID scheme and self-force 

To illustrate the entire scheme, consider one-dimensional (ID) motion of a 
single particle in the given steady-state fields depending on x-coordinate only. 
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Assume zero transverse components of the electric field as well as a zero longitu- 
dinal component of the magnetic field. Eqs.( |20|26l ) in ID read: 



where the subscript ± denotes the tranverse components of the vector potential 
and magneitc field and n^. is the vector along x-coordinate. In the limit of an 
infinitely small timestep and in negecting for a while the particle charge effect on 
the fields, we find that three integrals fully describe the mactroparticle motion: 

p(^) _|_ SE\i^)(_^Xp) = const, S + qp^{xp) = const, 

where, again, the continuous potentials at the particle position are expressed in 
terms of cell-centered grid function values. Within this approximation, however, 
the motion description is perfectly accurate, as long as the conservation laws are 
fulfilled which fully control the particle motion. 



However, as noticed in 11111 . 111311 . the energy conservaion property of the nu- 
merical scheme results in the side effect in a form of a self-force. The conserved 
energy includes the contribution from the particle's own field energy, which is 
weakly dependent on the particle location relative to the grid nodes. As a result 
the weak force appears tending to pull the particle from cell centers towards faces. 

Now one can calculate the self-force effect on the ID macroparticle motion. 
The current through x-faces, (^Filxp'^^) — Fj(Xp)) , can be integrated over 

time, giving the electric field: E^^^ = ^'^'^^^^ (Fi{xp) — Now we apply the 



weight coefficients as in Eg. (1201) to obtain the self-force acting on the macroparti- 
cle in the form of the gradient of the particle's own field energy: 



qpY.Mxp)E^^ 



Ax dxp 



With this accounting, the energy conservation law provides the distorted result: 

£ + qp (^{xp) + $'^°"'")(xp)) = const. Upon calculating the transverse magnetic 
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field and vector potential induced by the particle slow motion, one can find that 
the transverse particle momentum is distorted accordingly, so that the ratio, which 
is the transverse velocity, is not modified with the particle's own field. So, in ID 
non-relativistic motion the self-force does not affect the transverse velocity. 

To evaluate the self-force effect on the longitudinal velocity, note that the di- 
mensional factor in the formula for qp^^""^"^' ( related to temperature equals 
A^x / (Nrjj), where r/j is the Debye length and is the number of macroparticles 
per cell. Usually (but not in the numerical tests discussed below) this multiplier 
is small, which reduces the effect of the self-force comparing to the macropar- 
ticle interaction with the regular electromagnetic field and with the neighboring 



particles. The dimensionless factor was evaluated in Ullll for the lowest-order 
form-factor and the discontinuous character of the self-force was noticed for this 
case, resulting in potential problems. However, for the choice of a higher order 
form-factor (which is called "form-factor of the order / = 1^" below), the field 
energy is not only a smooth function of Xp, but its variance is very small: 

E = l+li^P-jnj+i- xp)\ j<xp<j+i. 

i 

The variance from minimum to maximum is 1/64. We see that the self-force is 
rather small and purely potential and its effect was not revealed in the numerical 
tests discussed below. 

4.5. Alternating-order form-factor 

Eqs. (ll4l20l26l) may be interpreted in terms of the alternating-order form- 
factor, Lp{x, Xp)(p{y, yp)Lp{z, Zp), with the value of this function at the grid point, 
X, y, z, giving the interpolation weight for the grid function defined at this point. 
Although the form-factor is a continuous function, only semi-integer and inte- 
gers X, y, z matter, to which the grid functions are assigned. Thus, we define 
^{x,Xp) = fi{xp) = f^^\i,Xp) at integer X, hut ip{x,Xp) = AFi+i/2(xp) = 
_l_ 1/2, Xp) at semi-integer x, i.e. we alternate the form-factor order. As 
long as the alternating-order form-factor combines the values of functions, f^'-^ 
and f^''~^^\ we would characterize it by the semi-integer order, = / + 1/2 

Finally, we provide the values of the alternating-order form-factor for = l|: 

z = int(xp), d = Xp-i, fi.,i+i{xp) = {1 - d; d), 

AFi_i/2:i+3/2{Xp) = 
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Figure 1 : Alternating order form-factors (shape functions) of the order 1 i (left panel) and 2 i (right 
panel), for different values of Xp, which are parameterized via d. The point values are presented 
for integer values of .t = i (marked with the vertical lines) and for semi-integer values of x. 

and for /a = 2|, which really ensures high accuracy and good energy conserva- 
tion: 



i = mi[Xp + -), d = Xp + --i, = I '^l-y2~''Y 



AFi_3/2:i+3/2{Xp) 



5. Numerical tests 



[1-df {2-df 2{l -df {l + df 2d^ d^' 



6 ' 6 3 ' 6 3 ' 6 



The tests we present combine some features of the tests for ChCPIC as de- 
veloped in ii. The tests are highly demanding: the mesh size is as large as 
the Debye length, the time step is as large as a half of Ax/c, the macroparticle 
density is as small as two per cell in a three-dimensional case, thus increasing 
the macroparticle charge. Increased effects of particle-particle and particle-grid 
interactions allow for producing a noticeable numerical error within a short test 
runs. 

First, we study the noise in the plasma with comparatively high level of ther- 
mal fluctuations. We use a three-dimensional 64*64*64 rectangular grid. The 
average number of electron particles per grid cell is = 2. The electron thermal 



speed is: vxe = yTe/me = 0.05c, where Tg is the electron temperature. At the 
initial time instant 2 *64^ electron particles are randomly distributed over the com- 
putational domain, the averaged plasma frequency being equal to Upe = 1 . The 
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Figure 2: Evolution of the total energy normalized per the initial total energy. Blue curves demon- 
strate perfect energy conservation as achieved with the alternating-order form-factor of the order 
la ~ 1 5 • The pink curves show much worse result obtained with higher (but not alternating) order 
I = 2. The left panel is for the Maxwellian plasma simulated through two plasma wave periods. 
The right panel is for two-stream instability followed till ojpet = 200 (8000 iterations). 



equal number of immovable ions are put to the same locations as the electrons, 
so that the electromagnetic field is zero initially. The mesh size is Ax = Ay = 
Az = vxe/^pe- With the time step, cOpeAt = 0.025, the plasma evolution has been 
simulated through two plasma wave periods, ncOpeAt ^ Air, n = 503. Fortran 
90/95 code is compiled and run with the double precision. 

In Fig. 2 (left panel) we present the evolution of the total energy normalized per 
the initial total energy. The blue curve demonstrates an almost perfect energy con- 
servation as achieved with the alternating-order form-factor of the order = l|. 
The pink curve shows, for comparison, the result obtained with the "conventional" 
way to interpolate the fields in the particle location, as used, for example in ^ 
as well as in older works, including the TRISTAN code Isj]. Within this approach 
the cell-centered form-factor provides the interpolation weight, which is applied 
to the electric field vector as averaged over the faces of the given cell, which re- 
sults in the following expression for a; 

</+i/2,fe+i/2(xp) = 1^ AFj+i/2(|/p)AFfc+i/2(^p), ... (31) 



and the analogous principle is applied to construct /^-coefficients, with Eg. (1181) 
being used to interpolate the electric currents. With this approach higher order 
/ = 2 is used to interpolate both the charge density and the electric and magnetic 
field. 

We clearly see that the increase in the interpolation order for the fields, com- 
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paring to the alternating-order form-factor results in the degrade of an accuracy, 
as long as the energy defect becomes two orders of magnitude higher. This numer- 
ical result entirely agrees with our theoretical claims that within ChCPIC scheme 
the energy conservation degrades unless the alternating-order form-factor is used. 

Note also that the worse quality of the results with the uniform order / = 2 
form-factor is accompanied by the degrade in efficiency. The stencil for / = 2 
form-factor is wider than that for la = l|, therefore, more search/algebra op- 
erations are involved in this case. Thus, the efficiency drops from ^ 2.1 ■ 10^ 
particles advanced per second of CPU time per processor, for = l| case, down 
to^ 1.5 ■ 10^ particles/s/processor, for / = 2 case. 

In the second test taken from ^ the two-stream instability is studied. The ini- 
tial distribution is modified as follows. First, the magnetic field is applied along 
the direction of x-axis, such that the cyclotron frequency of electrons in this field 
equals one: lubc = QeB/{mec) = 1 = cupe. To set the unstable distribution func- 
tion with two streams, for 50% of the electron particles the directed stream veloc- 
ity, V = 2vTe = O.lc, along x-direction is added to the thermal random velocity. 
For the other half of the electron particles the negative of the first stream velocity, 
—V = —2vTe = —O.lc, is added. The instability evolution is traced through the 
time interval, < copet < 200. In the right panel of Fig. 2, the check of the total 
energy conservation, again, demonstrates the advantage of the alternating-order 
form-factor (herewith the contribution from the initial magnetic field is excluded 
from the total energy integral). 

In Fig. 3 we present two test results obtained with the increased particle charge, 
so high that = 1 particle per cell simulates the plasma of unity plasma fre- 
quency cOpe = 1, i.e. the charge is two times higher compared to the previous 
tests. Left panel presents the Langmuir oscillations. The initial distribution is 
simular to that for studying the two-stream instability, but the number of particles 
(of two times higher charge) is by a factor of 0.5 less. The unidirectional stream 
velocity, v = 2vTe = O.lc, is added to all electron particles. The left panel of Fig. 3 
shows the total energy of electrons (upper curve) and the electric field energy, as 
functions of cOpet. The test parameters (such as the low number of particles, re- 
sulting in comparatively strong electron-ion correlations) are intentionally chosen 
to make the errors noticeable, so we can observe gradual attenuation of the oscil- 
lations, as the result of electron-ion interactions. However, the defect in the total 
energy is as small as ~ 10"^, thanks to using the alternating order from-factor. 
This defect is negligible compared to the changes in the partial energy integrals. 
In the right panel, we show the trajectories of particles gyrating in the magnetic 
field. The magnetic field is increased ten times: cobc = 10, and it is applied along 
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Figure 3: Test results with increased particle charge. Left panel: particle energy and longitudinal 
electric field energy for the Langmuir oscillation. Right panel: electron and positron trajectories 
in the magnetic field, black lines show the grid (cell boundaries). 

2;-axis. Two particles are located initially in the center of the computational do- 
main, with oppositely directed initial momenta, w = ±0.5c, the particle charges 
have opposite signs (electron and positron particles). The particle trajectories are 
shown in the right panel of Fig. 3 for the time interval, < to Bet < 400 with the 
timestep, cose^t = 0.25. One can observe a small gradual decrease in the parti- 
cle kinetic energy. It reveals itself in the curve "thickness", with its inner radius, 
which is the Larmor radius at the final particle energy, being somewhat less than 
the outer radius, which is the Larmor radius at the initial particle energy. 

Our comparison with the field interpolation scheme presented in Eq. (|3TI) does 
not aim to criticize [Isl, 0], moreover that Eq.(l3TI) in TRISTAN is used in a com- 
bination with low order / = 1 form-factor for the charge density, so that reducing 
this order following Eq.(|20l). would be impractical and still would not allow us 
to reach the energy conservation. However, with the higher order form-factor, as 
we see, the development of the interpolation scheme should go along the direction 
different from that presented in Eq.(l3T1). 

6. Conclusion 

The pulsed electric field in a focus of a high-field laser may cause a charge 
separation even in the target of a solid-state density. To simulate the laser-plasma 
interaction in that strong fields, the ChCPIC scheme could be the best choice (see 
10]). However, these schemes are believed to be too noisy and the energy non- 
conservation, as we see, can potentially be a source for such noise. 
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We show that the ahemating-order form-factor is to mitigate this inherent flaw. 
I am grateful to Prof. T. Zh. Esirkepov, Prof. K. Powell and Dr. N. Naumova 
for discussions. 



Appendix. Full algorithm to advance the PIC solution. 

Assume the magnetic field and the particle momenta to be known att = n — 
1/2, as well as the electric field and particle positions att = n. Advance all the 
fields and particle through a time step. First, advance the magnetic field through a 
half time step: 

rjW n _ TDix) n-l/2_ 
-°i+l/2,i,fc - ^i+l/2,j,k 

2 V i+l/2,i+l/2,fc -^i+l/2,j-l/2,fc/'"r 2 «+l/2,i,A:+l/2 -^i+l/2,j,k-l/2) ^ ■ ■ ■ 

Alternatively, the vector potential may be updated through a half time step. 

Then the particle motion is updated. For each particle, first, the fields at the 
particle position should be interpolated using Eq.(l20l) and either Eq.(l26l) with the 
updated magnetic field or Eq.(l24l) with the update vector potential. To do this, the 
alternating-order form-factor should be calculated for the particle position, x^. 
Then, the particle momentum is advanced through a half time step, accounting for 
the effect from the electric field only: 

w; = w;-V2 + e(x;), e(x;) = |^E(x;), 

Then, the contribution from the magnetic force is added: 

w," = < + [w; X b(x^)], b(x^") = ^^==B{^;). 

The momentum is advanced through the full time step then: 

^1/9 , . 2[w^xb(x")] 

n,+l/2 = I eCx'^) I \^ 

the correction in the last term is chosen in such a manner that the magnetic field 
does not affect the particle energy, so that 

(w^^/' - e(x;))2 = (w^"-V2 + e(x^-))^ (33) 
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Finally, the particle current should be calculated. First, save the form-factor 
at the cell centers, AFi+i/2(Xp), AFj+i/2(|/p), AFk+i/2{Zp). Then calculate the 
particle velocity and new particle position: 

^n+l/2 ^"+1/2 

^ ^ x;+i = x; + ^^ diag(c,, c„c,) (34) 



Find the new form-factors, AF,+i/2(x^+^), AF,+i/2(y^+^), AFfc+i/2(^^+^), then 
calculate AF.^^^^g' ^^^^ calculating the partial sums of 

^-^i+i/2' ^-^^-+1/2' ■^-^fc+i/2' fi'^^ '^he currents, ^, using Eq.(fT8l). End update for 
this particle and proceed to the next one. 

The magnetic field then should be advanced through another half time in the 
way as presented in Eq.(l32l). Finally the electric field should be updated, with the 
electric current density, which sums up the contributions as in Eq.(fT8l) from all 
particles: 

rp{x) n+1 _ rp{x) n _ ax, j{x) n+1/2 

-^iJ+l/2,fc+l/2 ~ -^ij+l/2,fc+l/2 ^''^''"'j,j+l/2,A,'+l/2"'~ 

, ^ CR^'') _ n+1/2. _ ,T^{y) n+1/2 _ n+1/2. 

^y\^i,j+l,k+l/2 ^i,j,k+l/2 J ^zyJ->ij+i/2,k+l ^i,j+l/2,k )^ ■■■ W^'-' 

End update for this time step and proceed to the next one. 
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